function p = mda(X, lab, d)
% Perform MDA on the data set X with the class labels in lab
% [p] = mda(X, lab)
% lab is a column vector containing the class labels, which are 
% assumed to be 1, 2, ... K. K is the number of class.
% Each row of X corresponds to a pattern. So, for the entire 
% Iris data set, X will of size 150 by 4.
% lab in this case should be of size 150 by 1 with entries from 1, 2, 3.
% d is the dimension to be reduced.
% Note that d < K. (Why?)
%
% For output,
%   p is the direction for projection; each column is one direction
% So, suppose you have another set of pattern stored in the matrix Y,
%   where each row in Y corresponds to a pattern.
% You should reduce the dimensionality of Y by:
% 
%   yy = ( Y * p )';
%
% Written by Martin law lawhiu@cse.msu.edu
% Jan 6 2005

K  = max(lab);
D  = size(X,2);

mean_all = sum(X,1) / size(X,1);
m  = zeros(K, D);
n  = zeros(K, 1);
Sw = zeros(D);
Sb = zeros(D);
for ii=1:K,
  current = X(lab == ii, :);
  n(ii)   = size(current,1);
  m(ii,:) = sum(current,1) / n(ii);
  current = current - repmat(m(ii,:), n(ii), 1);
  Sw = Sw + current' * current; 
  temp = m(ii,:) - mean_all;
  Sb = Sb + n(ii) * temp' * temp;
end  

opt.disp = 0;
[p, D] = eigs(Sb, Sw, d, 'LM', opt);

for ii=1:d,
  p(:,ii) = p(:,ii) / sqrt(sum(p(:,ii).^2));
end

return;
